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We present a new algorithm for quantum Monte Carlo simulation based on global updating 
with loops. While various theoretical predictions are confirmed in one dimension, we find, for 
S = 1 systems on a square lattice with an antiferromagnetic biquadratic interaction, that the 
intermediate phase between the antiferromagnetic and the ferromagnetic phases is disordered and 
that the two phase transitions are both of the first order in contrast to the one-dimensional case. 
It is strongly suggested that the transition points coincide with those at which the algorithm 
changes qualitatively. 
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Most theoretical studies on the model described by the 
Hamiltonian 



Jq < 0. Provided that Jq > 
available for the model when 



0, the new algorithm is 



-H = J2(J L {S i -S j ) + J Q (S i -S j f) (1) (2S 2 -25+1) < J L /J Q or J L /J Q < -2S(S-l), (2) 



have been lead mainly by purely theoretical motiva- 
tions. There are many special points in the phase di- 
agram where one can obtain some rigor 
the model with S = 1 in one dimensii 
In particular, it is generally believed that at each phase 
transition point, the model is integrable and a Bcthc 
ansatz solution exists. Recent experiments on a quasi 
one-dimensional system, LiVGe2 06,c!P however, pro- 
vided us with a new motivation for studying this model. 
The experimental results suggested the presence of a bi- 
quadratic interaction of considerable magnitude. The 
succeeding numerical workliiP indicated that the magni- 
tude was not sufficiently large to take the system out 
of the Haldane phase. However, it also confirmed that 
the biquadratic interaction considerably affects the na- 
ture of the system. Since there is no reason to expect a 
non-negligible contribution of the biquadratic interaction 
only in the one-dimensional case, it seems reasonable to 
start studying the system in higher dimensions. 

Unlike the one-dimensional case, our understanding of 
systems in higher dimensions is limited, due to a lack of 
exact solutions and powerful field-theoretical methods. 
In this paper we investigate the model with a Monte 
Carlo method using a newj-algorithm based on the con- 
cept of the loop algorithm.t3 It is found that the phase 
transition points for the S = 1 case in two dimensions 
coincide with (or are located very close to) the points 
where the types of graphs in the loop algorithm change 
as shown below. 

First, it should be noted that, similar to the con- 
ventional world-line algorithm, the new algorithm also 
encounters the negative sign difficulty with the ferro- 
magnetic biquadratic interaction, that is, in the case of 



for general spin S and in arbitrary dimensions. Outside 
this region, we encounter the negative sign difficulty. In 
what follows we do not discuss algorithms in such cases, 
since they would not be of much use. In contrast to the 
case of frustrated spin models, the difficulty remains even 
if there are no closed (frustrated) cycles in the system. 
We encounter serious difficulty even in one dimension. It 
is found, however, that the case of S = 1 is exceptional. 
As we see below, we are exempt from negative signs as 
long as Jq > holds. 

The method we use here is the world-line Monte Carlo 
method. Recently, a number of new techniques have been 
developed for reducing the long computational correla- 
tion times that are problematic in many applications. 
Loon algorithms were found to be particularly success- 
fulOtj A loop algorithm consists of two procedures; 
one for constructing a graph (graph assignment) and the 
other for generating a configuration (loop or cluster flip) . 
In general, both procedures are carried out probabilisti- 
cally. The graph assignment is applied to each plaquette 
on nearest-neighbor sites and short (formally infinitesi- 
mal) intervals of imaginary time. In simple cases, we can 
easily compute the graph assignment probability once we 
express the Hamiltonian in the following form, 



H 



local 



E 

Ger 



o(G)A(G) (a(G) > 0), 



(3) 



where rii oca \ is the local Hamiltonian and the right-hand 
side is a sum of operators over a certain set of graphs 
T. The symbol A(G) denotes an operator whose matrix 
elements are one when the graph G matches the initial 
and the final spin states, and zero otherwise. The proba- 
bility for assigning a graph G to an interval of the length 
At is simply a(G)Ar if the matrix element of A(G) cor- 
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Fig. 1. The four types of graphs. The vertical direction corre- 
sponds to the imaginary time. Each pair of grey vertical repre- 
sents two Pauli spins that constitute one of the original spins, in 
the case of S = 1. 



responding to its current state is one. When no graph 
is assigned for a plaquette, we assign the identity graph 
which is simply a pair of straight lines parallel to the 
time axis. After all graph assignments, each loop which 
consists of graph-linked spins is independently flipped. 

In the case of 5 > 1/2, we express each spin operator 
as a sum of 25 Pauli operators,!! 2 }^ 

^ 2S 

Si => - 22 <7i,n, 

rendering the Hilbert space greater than the original 
one. In order to eliminate contributions from unnec- 
essary states, we insert a symmetrization operator P at 
t = (3. This operator projects the entire phase space into 
the original space in which E^i^/^ 2 = S(S + 1) 
for all i. In addition, when we know that the negative 
signs always appear in pairs that cancel each other, we 
can simply replace negative off-diagonal elements of the 
Hamiltonian by their absolute values. Thus, in general, 
we can replace eq. (^|) by 



[| — 7Ylocal|]s 



]T a(G)A(G) ] 
,Gf=r / 



(4) 



where a(G) are positive numbers and \Q\ denotes an 
operator whose off-diagonal matrix elements are equal 
to the absolute values of corresponding elements of Q. 
The symbol [• ■ -] s denotes the symmetrization, that is, 

The following identities are useful for obtaining an ex- 
pression of the form of eq. (||) in the present case. 

[S t - Si], = 5 2 (-[l] s + 2A sc ), (5) 



[{Si ■ S/f} s = 5 2 ((5 2 + l)[l] s - 2(2S 2 - 2S - 
+ (25-l) 2 A DC ). 



1)A SC 
(6) 



The operator Asc is defined to be [A(Gsc)]« with anv 
single cross graph Gsc (Fig- 1)- Clearly, the definition 
does not depend on the choice of Gsc because of the 
symmetrization. The other operator Adc> is defined 
similarly with a double cross graph Gdc- Then, it is 
easy to obtain the following expression 

[|-Wlocal|] s = [\MSi ■ Sj) + J Q (Si ■ S 3 ) 2 \] s 



2\K L -K Q 



2S 2 -25+1 
(25-1) 2 



A. 



sc 



^qA dc 



(7) 



where K L = S 2 J L and K Q = S 2 {2S - 1) 2 Jq. Unim- 
portant constant terms have been omitted. If all the 
coefficients are non-negative, this equation already has 
the form of eq. @. Therefore, if Kq > and K l /Kq > 
(25 2 -25+l)/(25-l) 2 , eq. (7) defines a valid algorithm. 
This algorithm is termed the ferromagnetic algorithm in 
the following. 

If Jl < and \Jl\ is sufficiently large, we do not have 
negative sign configurations. We can see this easily using 
the unitary transformation of spin operators on one of 
two sublattices: S x -> -S x , S y -> -S y and S z -> S z . 
It is equivalent to replacing Si ■ Sj by Si a Sj = —SfS x — 



SfS v 



, , SfSj which yields the following equation. 



— Wlocal|]s = [\JL(Si O Sj 

2S(S - 1 



-Kr 



(2s-iy 



J Q {Si o Sj)% 
K Q ) A SH + ^qAdh 



,(8) 



where Ash and Adh are the operators corresponding 
to the single horizontal graph and the double horizontal 
graph (Fig. 1), respectively, analogous to Asc and Adc 
mentioned above. The coefficients are all non-negative if 
K Q > and K l /Kq < -25(5 - l)/(25 - l) 2 . In this 
case, the above expression is of the form of eq. (Q) pro- 
viding us with a valid algorithm. We call this algorithm 
antiferromagnetic. 

On the other hand, in the intermediate region —25(5— 
l)/(25-l) 2 < K l /Kq < (25 2 -25+l)/(25-l) 2 , we 
cannot avoid the negative matrix elements in general. 
However, in the case of 5 = 1, we see that the negative 
signs always cancel each other to give a positive global 
configuration. This is because only the matrix elements 
between a state with two holes (the site with 5 Z = 0) and 
a state with antiparallel spins are negative, and all the 
others are positive. Since the number of holes changes 
only in transitions corresponding to this type of matrix 
element, the number of pair creations of holes must be 
equal to that of pair annihilation. In order to construct 
a loop algorithm in this case, we note that 



-H 



local |Js — 



(Kq — Kl) Adh + ^l^dc • (9) 



This gives a set of probabilities required for the algorithm 
called intermediate. 

It is worth noting that in the case of 5 = 1, both of 
the algorithmic transition points correspond to special 
points in the physics of the model. Using standard pa- 
rameterization, we define 6 and J so that Kl = — J cos 9 
and Kq — — J sin 9. Then, the algorithmic transition 
from the ferromagnetic algorithm to the intermediate al- 
gorithm occurs at 9 — —n/2, whereas the transition from 
the intermediate algorithm to the antiferromagnetic one 
occurs at 6 = —Sir/A. The latter (9 — — 3tt/4) is known 
to be the real transition point in one dimension at which 
the non-magnetic (dimer) ground state switches to the 
ferromagnetic ground state. The former (9 = —n/2) 
corresponds to the point at which a Bcthe ansatz solu- 
tionO'EfP is available. At this particular point, the algo- 
rithm consists of only one type of graph, i.e., the double 
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horizontal graphs. In terms of graphs, the system is very 
similar to two independent anti-ferromagnetic Heisen- 
berg models with S — 1/2. The only difference is that 
the two systems are coupled with each other through 
the symmetrization operator at r = (3. The situation at 
9 = —3tt/4 is similar. In graphical terms, at 9 — — 37r/4, 
the system is two 5=1/2 Hciscnbcrg ferromagnets cou- 
pled to each other only at r = (3. We suspect that there 
is a relationship between these facts and the integrability 
of the system at these points. 

It should also be noted that the locations of the algo- 
rithmic transition points do not depend on the dimen- 
sionality or the lattice structure since they only depend 
on properties of the two-point Hamiltonian. As we will 
see below, the present numerical results strongly suggest 
that the algorithmic transition points arc special not only 
in terms of the algorithm but also in real physics even in 
higher dimensions. 

To evaluate the validity of the method, we have veri- 
fied various predictions based on exact solutions in one 
dimension. We first locate the phase transition points. 
It is easy to do so for the dimer-ferromagnetic phase 
transition at 9 — — 37r/4 since there is a very clear dis- 
continuity there in the first derivative of the energy with 
respect to 9, indicating that the ground state switches 
from the dimer state to the ferromagnetic state. In ad- 
dition, there is a clear discontinuity in the magnetization 
as a function of 9 at 9 = — 37r/4. 

On the other hand, it is difficult to see a singularity in 
energy as a function of 9 around 9 = — ir/4. We compute 
the Binder cumulant ratio 

for the staggered magnetization and the dimer order pa- 
rameter. Here, Q is the operator of an observable as 
given above and {{Q n )) denotes the temporal average as 
well as the thermal average, i.e., 

((Q")) = JL f dn • • • dT n T(Q( Tl ) ■ ■ ■ Q(r„)), (11) 
P Jo 

where T denotes the time-ordered product and the single 
bracket denotes the thermal average. We have seen that, 
with a fixed value of (3/L = 1/2 for example, gM s curves 
(as a function of 9) for various L intersect one another 
at 9 = — it /A. In addition, when plotted as a function 
of (3 at 9 = — 7r/4, a curve for a fixed system size has a 
peak around (3 ~ /3 pca k- The location of the peak /3 poa k is 
proportional to L and the peak height is independent of 
the size, indicating that 9 = — 7r/4 is critical and z = 1, 
as expected. We also compute the correction terms in 
E(L,T = 0) expanded with respected to 1/L 2 and in 
C(L = oo, T) expanded with respect to T. From the 
coefficients of these we have estimated the central charge 
and the velocity. They are found to be consistent with 
the theoretical predictions c = 1 and v = n/3. 

Now we turn to the two-dimensional case. We have 
performed simulations for up to L — 64 and (3 = 32 at 
various values of 9. First, similar to the one-dimensional 
case, we can easily locate the phase transition to the fer- 
romagnetic phase in the energy curve as a function of 
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Fig. 2. Total energy in the case of S = 1 in two dimensions. 
The inset is a magnified view of the vicinity of the algorithmic 
transition point 9 = —n/2. The data points are for sufficiently 
large L and {3 so that they have no visible dependence on L and 
0- 
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Fig. 3. Staggered magnetization per spin at the algorithmic tran- 
sition point, 9 = — 7r/2. 



9 (Fig. 2). It is easy to prove that the ground state 
is purely ferromagnetic in the region 9 < — 37r/4. The 
present result suggests that as soon as the simple proof 
becomes invalid, the ground state switches from the fer- 
romagnetic state to another state with different symme- 
try, in the same way as the one-dimensional case. We 
can also observe the discontinuity in the magnetization 
curve. Thus, we see that one of the phase transitions co- 
incides with one of the two algorithmic transition points. 

For the other transition point, we compute the stag- 
gered magnetization. Since the system has a finite stag- 
gered magnetization at 9 — 0, we expect at least a transi- 
tion from the antiferromagnetic phase to another phase. 

We observe that, for each system size, the staggered 
magnetization increases initially and saturates at its 
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Fig. 4. Staggered magnetization per spin in the thermodynamic 
limit. 



zero-temperature value as the temperature is lowered. 
As an example, the staggered magnetization for various 
system sizes and temperatures is plotted at 9 — — 7r/2, 
in Fig. 3. To the zero-temperature values for various 
system sizes we fit the function 



M s (£, 9) /L 2 = M s (oo, 9) /L 2 + aL 



-d/2 



in order to obtain the values in the thermodynamic limit. 
In this way we estimate the staggered magnetization per 
spin at the algorithmic transition point, i.e., 9 = —tt/2 
as 



M s {6 = -n/2)/L 2 



0.14. 



(12) 



This value is small but finite, whereas for 9 = — 0.505-7T, 
we find the staggered magnetization is vanishing. There- 
fore, one of the following two possibilities is correct. One 
possibility is that the antiferromagnetic-paramagnetic 
transition is of the second order and very close to but 
slightly below —tt/2, and the other is that the transition 
is of the first order and occurs somewhere in the interval 
(-0.5057T, -0.5007T] including 9 = -tt/2. 



To summarize, we have presented a new loop algorithm 
for the Heisenberg model with the biquadratic term. In 
the case of S = 1 it covers the entire region of 9 < 0, 
while it covers a part of it for S > 1. Outside these re- 
gions, not only the loop algorithm but also the general 
world-line Monte Carlo simulation suffers from the neg- 
ative sign difficulty. The algorithmic transition points 
have been found to be special from the physical point 
of view. In the two-dimensional case, in particular, they 
correspond to the real transition points. In contrast to 
the one-dimensional case, the intermediate phase is dis- 
ordered. A natural interpretation of the present results 
also suggests that the transition into this phase from the 
antiferromagnetic phase is of the first order. The ex- 
tension of the algorithm to the cases where the simple 
methods encounter the negative signs is an open prob- 
lem. 
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In Fig. 4, we plot M s / L 2 in the thermodynamic limit 
(extrapolated first to (3 = oo and then to L — oo) as a 
function of 9. It seems that, as one approaches this point 
from above, M s /L 2 converges to a finite value different 
from the one in eq. (^)- We estimate this limiting value 
as 

lim M s (9)/L 2 ~ 0.24. (13) 

This suggests that the latter of the two possibilities is 
correct and that the transition is exactly at the algorith- 
mic transition point. 

Regarding the intermediate phase — 37r/4 < 9 < —tt/2, 
we can see no evidence of any kind of order. At least we 
see that there is no long-range correlation in the energy, 
the magnetization, the staggered magnetization, or the 
(tt, 0) or (n, 7r) structure factor of the dimer order pa- 
rameter. 



